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SECTION  1 


INTRODUCTION 


In  this  report  we  shall  describe  a  method  for  calculating  the  propagation 
of  elastic  waves.  The  method,  the  phase-screen  method,  has  some  relation  to  WKBJ 
and  other  high  frequency  methods;  the  details  are  quite  different.  The  wave  equations 
of  elastodynamics  are  solved  exactly,  the  approximation  being  made  on  the  medium 
through  which  the  waves  propagate  rather  than  on  the  solutions  to  the  dynamical 
equations  themselves.  While  the  method  can  be  used  to  calculate  the  propagation 
through  a  complex  medium  with  a  known  spatial  distribution  of  elastic  parameters 
(an  example  is  given  below),  probably  the  most  common  use  of  the  method  in  previous 
work  has  been  to  study  the  propagation  of  waves  in  stochastic  media  where  only 
statistical  properties  of  the  media  are  known;  it  is  this  last  application  we  have  in 
mind  here. 

The  method  is  not  new.  It  has  been  used  for  many  years  by  astronomers 
to  estimate  the  effects  of  the  atmosphere  on  the  propagation  of  starlight,  Ratcliffe 
[1956],  by  engineers  to  study  the  effects  of  structure  on  the  propagation  of  signals 
from  communication  satellites,  Knepp  [1983],  by  acoustical  engineers  to  study  the 
propagation  of  sound  in  the  ocean,  Flattd  [1983]  and  Martin  and  Flattd  [1988],  and 
for  other  purposes.  All  of  the  past  studies  of  which  we  arc  aware  have  considered 
only  the  case  of  scalar  waves.  Here  we  shall  formulate  the  problem  for  the  case  of 
vector  elastic  waves;  indeed  it  is  the  interconversion  of  transverse  and  longitudinal 
components  that  provides  a  major  focus  of  the  study. 

Consider  a  simple  plane  wave  ( k  -  tu/c),  where  a  harmonic  time  dependence 
of  frequency  w  is  assumed  propagating  in  a  medium  with  sound  speed  e. 

.w 

i— a  .  . 

^  *  c  •  (1) 

If  we  use  that  solution  as  an  approximation  in  a  medium  with  sound  speed  e\  we  find 
that  at  any  fixed  a  the  relation  between  the  approximate  solution,  4,  and  the  exact 
solution,  $ ,  is  just  a  phase 


In  this  simplest  example  the  effect  of  a  velocity  anomaly  is  just  to  advance  or  retard 
the  phase  of  the  wave;  that  is  the  central  approximation  in  the  phase-screen  method: 
we  assume  that  the  effect  of  velocity  anomalies  is  to  retard  or  advance  the  phase  of 
the  wave. 

Now  consider  the  situation  depicted  in  figure  1.  A  wave  is  propagating  from 
z  =  0  to  z  =  As  through  a  medium  whose  mean  sound  speed  is  e  with  velocity  anoma¬ 
lies  6c(x,y,z)  present.  Our  first  approximation  is  that  the  effect  of  the  anomalies  is 
to  multiply  the  solution  at  As  by  an  x-y-dependent  phase  factor  given  by 


exp  (*A(x,y)J  =  exp 


.to 

-i— As  -I-  tto 
c 


dz' 

c  +  6c(x,y,z') 


(3) 


Thus  the  problem  shown  in  figure  1  has  been  replaced  by  that  shown  in  figure  2  where 
the  propagation  from  0  to  z  -  As.  is  through  a  medium  with  a  constant  sound  speed 
c.  The  result  is  then  multiplied  by  a  phase  factor  given  by  (3)  to  obtain  the  wave  at 
s  =  As+. 


To  this  point  we  have  considered  the  case  where  at  z  =  0  we  had  a  simple 
plane  wave  but  for  more  complex  cases  where  the  solution  at  z  =  0  is  a  function  of 
x  and  y  we  will  use  the  same  procedure  of  figure  2:  propagate  through  a  uniform 
medium  to  s  =  As.  then  correct  the  phase  at  s  -  Asf .  With  that  generalisation 
we  have  the  possibility  of  repeating  the  procedure,  inserting  phase  screens  at  other 
larger  values  of  s.  That  possibility  brings  up  the  other  approximation  we  must  make: 
we  must  make  the  problem  parabolic;  that  is  to  say,  we  wish  to  have  a  formulation 
such  that  the  value  of  the  wave  function  at  some  fixed  s  determines  the  value  at 
larger  s’s.  Since  the  wave  equation  is  second  order  in  each  spatial  variable  that  is 
not  the  usual  case;  usually  an  initial  value  problem  would  require  specification  of  two 
functions  (e.g.  d  and  d(d)  at  a  fixed  z  in  order  to  determine  the  solution.  As  an 
alternative  we  can  specify  the  value  of  the  function  on  two  planes  corresponding  to 
different  values  of  z-  What  we  shall  do  is  a  variant  of  this  last  possibility:  we  shall 
use  the  value  of  the  function  at  z  ~  0  and  impose  outgoing  boundary  conditions  at 
a  -*  oo.  That  formulates  a  well-posed  boundary  value  problem  and  allows  us  to  solve 
for  the  function  at  As* .  A  standard  way  to  solve  such  a  problem  is  to  decompose 
the  value  of  the  function  at  s  0  into  a  set  of  plane  waves  whose  propagation  vector 
points  into  the  right-hand  cone.  Thus  we  make  the  problem  parabolic  by  a  procedure 
which  can  be  stated  one  of  two  equivalent  ways:  we  use  the  value  of  the  function  at 
iaO  plus  a  boundary  condition  at  s  -*  oo;  or,  we  use  the  value  of  the  function  at 
*  ss  0  plus  the  assumption  of  forward  propagating  wave*. 
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INTEGRATION  LINES 


Figure  1.  Propagation  of  the  wave  through  a  segment  of  inhomogeneous 
medium  alters  the  phase  of  the  wave  as  a  function  of  position 
due  to  the  velocity  anomalies  in  the  medium. 


Figure  3.  The  phase- screen  method  replaces  the  inhomogeneous  segment 
with  a  uniform  segment  and  a  pair  of  phase  screens.  The 
integrated  phase  defect  of  the  segment  is  prelected  onto  the 
second  screen. 
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The  approximations  just  described  allow  us  to  take  an  initial  value  for  the 
wave  function  and  propagate  the  wave  through  an  arbitrary  number  of  phase-screens 
to  find  its  value  at  any  2.  Furthermore  the  extension  to  vector  elastic  waves  is  obvious: 
given  the  displacements  at  2  -  0  we  can  uniquely  decompose  them  onto  the  set  of  S 
and  P  waves  which  have  propagation  vectors  pointing  into  the  right-hand  cone.  We  can 
thus  evaluate  the  displacements  at  z  -  A*_,  modify  them  at  2  =  A«+  with  separate 
phase  factors  using  the  appropriate  phase  velocity  for  each,  and  then  propagate  the 
wave,  according  to  the  uniform  vector  wave  equation,  to  the  next  phase  screen  where 
the  process  is  repeated.  Detailed  formulas  are  given  in  the  next  section. 


We  should  remark  here  about  the  relation  between  the  explanation  of  the 
phase-screen  method  we  have  just  presented  and  another,  slightly  different  discussion 
the  reader  may  have  seen  elsewhere.  Sometimes  the  starting  point  for  explanation  of 
the  phase-screen  method  is  to  write  the  field  as 


(4) 


where 


k 


w 

c 


(5) 


and  is  a  slowly  varying  function  of  s.  This  leads  to  the  replacement  of  the  wave 
equation 


* 


(e) 


with  the  “parabolic  wave  equation* 


V*  +  2ik, 


i h 


where 


(?) 


v\ «  %  +  a*.  (8) 
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The  problem  is  thus  made  parabolic  by  modifying  the  dynamical  equation  to  be  first 
order  in  a  so  specifying  4>  at  one  value  of  x  determines  the  solution. 

In  addition  to  this  modification  of  the  wave  equation,  the  phase  factors  that 
correct  for  the  velocity  anomalies  of  the  medium  are  generally  taken  to  be  independent 
of  k,  and  kv.  These  approximations  are  justified  for  either  scalar  or  vector  waves 
provided  k,  and  kv  are  in  fact  small  compared  to  w/c,  waere  £  genericalty  represents 
either  of  the  two  average  phase  velocities  for  the  vector  wave  case.  By  dimensional 
analysis,  kt  and  kv  are  inversely  proportional  to  the  long;,  scales  in  the  plane  normal 
to  the  r-axis.  Thus  if  these  lengths  are  large  enough,  .ae  approximation  is  valid. 
We  have  found,  based  on  the  comparison  of  the  phase-s  -een  method  with  an  exart 
solution,  which  we  describe  below,  that  using  the  f;,r  uniform  wave  equation,  and 
phase  factors  that  depend  on  J It,  and  lb,,  allows  us  to  ccurately  calculate  propagation 
for  slightly  smaller  structures.  Certainly  larger  structures  will  lead  to  better  results 
since  waves  propagating  at  large  oblique  angles  a.-e  apt  to  induce  backscaitering,  which 
is  not  accounted  for  in  the  formalism.  For  large  enough  structures,  thr.  parabolic 
approximation  is  recovered  by  expanding  cur  results  for  small 

A  slight  variant  of  the  method  has  been  used  before  in  seismology.  Haddon 
and  Huseby  { 1978)  used  what  amounts  to  a  phase-screen  calculation  combined  with  a 
stationary  phase  (ray)  approximation,  Mercier  [1961],  to  calculate  propagation  into  the 
NORSAIt  array.  The  combination  of  phase-screen  and  stationary  phase  approximation 
produces  the  thin  lease  approximation  of  geometrical  optics.  They  considered  only  P 
waves.  The  method  described  here  could  be  used  for  the  same  situations  including 
both  S  and  P  waves. 

Until  now  we  have  simply  generalised  -he  phass-screen  approximation  for 
scalar  waves  in  an  obvious  way  to  account  for  the  two  phase  velocities,  and  avoided 
making  the  parabolic  approximation  to  the  wave  equation.  Regardless  of  whether 
we  employ  the  parabolic  approximation  or  not,  this  leads  to  significant  consequences 
regarding  the  conservation  of  energy.  Certainly  if  a  wave  is  multiplied  by  an  overall 
phase,  energy  will  be  conserved.  However,  if  individual  contributions  to  the  total  wave 
are  multiplied  by  different  phases,  the  energy  of  the  subsequent  wave  will  not  be  equal 
to  the  energy  of  the  original  wave;  interference  terms  will  result.  Indeed  this  is  the 
case  if  we  naively  multiply  the  S  and  P  waves  by  different  phase  factors.  Although  the 
violation  of  energy  conservation  is  quite  small  at  each  phase  screen  using  parameters 
for  whkh  we  expect  the  method  to  be  valid,  the  cumulative  result  of  propagating  the 
wave  over  distances  of  say  1000  km  is  significant.  U  is  also  particularly  disconcerting 
that  energy  may  be  gained  as  well  as  lost. 
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To  remedy  this  problem  while  maintaining  the  spirit  of  the  method,  we  mod¬ 
ify  the  method  described  thus  far  in  the  following  way.  First,  we  find  the  recursion 
relations  that  express  the  Fourier  coefficients  of  the  wave  after  it  has  passed  through 
the  screen  in  terms  of  the  Fourier  coefficients  of  the  wave  before  the  screen  and  the 
phase  factors  which  modify  the  wave  at  the  screen.  The  recursion  relations  are  conve¬ 
niently  expressed  in  terms  of  a  matrix  multiplying  the  column  vector  of  old  coefficients 
to  produce  the  column  vector  of  new  coefficients.  These  are  vectors  and  matrices  in 
the  space  of  Fourier  modes,  not  spacetime.  Thus  if  V ,  V\  and  U  are  the  old  vector, 
the  new  vector,  and  the  recursion  matrix  respectively,  then 


v;  =  u^v.,  (9) 

where  we  are  employing  the  convention  that  repeated  indices  are  summed  over. 

In  principle  the  vector  space  we  are  dealing  with  is  infinite  dimensional, 
however,  for  practical  computational  purposes  it  must  be  truncated.  The  second  step 
of  the  algorithm  is  to  neglect  the  evanescent  modes,  keeping  only  the  finite  number 
of  modes  with  harmonic  ^-dependence.  This  is  reasonable  for  two  reasons.  First, 
we  are  interested  in  propagation  distances  for  which  the  evanescent  modes  have  been 
exponentially  damped.  And  second,  as  we  will  show  later,  the  evanescent  modes  do 
not  contribute  to  the  energy.  In  fact,  the  energy  flux  in  the  2-direction  averaged  over 
time  and  the  coordinates  perpindicular  to  the  2-axis  is  simply  the  norm  squared  of 
the  vector  of  Fourier  coefficients  corresponding  to  the  modes  that  propagate. 

The  condition  that  energy  is  conserved  across  the  phase  screen  is  that  U  be 
a  unitary  matrix,  t.  e.  it  must  preserve  the  length  of  the  vector.  This,  however,  is  not 
the  case  when  we  naively  multiply  the  S  and  P  waves  by  individual  phases,  although 
for  cases  of  interest  it  is  almost  unitary  in  the  sense  that  its  eigenvalues  are  almost  of 
unit  modulus,  and  the  commutator  of  U  with  its  adjoint  has  small  entries.  1  Thus 
the  third  step  of  the  algorithm  is  to  make  U  unitary.  The  procedure  used  to  modify 
the  matrix  will  be  described  later  in  full  detail. 

We  are  now  ready  to  move  on  to  the  details  of  the  method,  which  will  be 
organized  as  follows.  In  chapter  2,  we  present  the  general  theoretical  details  of  the 
3D  method.  In  particular,  the  general  recursion  relations  that  determine  the  Fourier 
coefficients  of  the  succeeding  wave  in  terms  of  the  Fourier  coefficients  of  the  previous 
wave  and  the  A -phases, ».  e.  the  phases  that  modify  the  wave  at  the  phase  screen,  are 

lThe  necessary  and  sufficient  conditions  for  a  matrix  to  be  unitary  is  that  it  commute  with  its 
adjoint  and  have  unit  modulus  eigenvalues. 


derived.  With  these  recursion  relations,  an  incoming  source  wave  may  be  iteratively 
propagated  through  an  arbitrary  number  of  phase  screens. 

Also  in  chapter  2,  we  derive  a  general  criterion  to  specify  the  number  and 
spacing  of  phase  screens  to  be  used.  The  basic  idea  behind  the  criterion  is  that  for  fixed 
observation  distance  the  phase-screen  solution  should  converge  to  a  fixed  answer  as  the 
number  of  phase  screens  are  increased.  This  is  analogous  to  the  calculus  of  integrals 
where  as  the  interval  is  broken  up  into  smaller  and  smaller  regions  the  convergent  limit 
is  obtained.  2  To  optimize  computer  time,  we  are  interested  in  finding  the  marginal 
number  of  phase  screens  needed  such  that  the  answer  converges  to  its  asymptotic  limit 
to  within  a  fixed  tolerance.  The  criterion  is  cast  in  the  form  of  a  Cauchy  condition 
for  sequences;  two  successive  terms  in  the  sequence  are  compared  to  determine  how 
well  the  sequence  is  converging. 

Finally  in  chapter  2,  the  time-averaged  energy  flux  is  computed  in  terms 
of  the  Fourier  coefficients  of  the  phase-screen  solution.  As  mentioned  previously  the 
displacement  decomposes  naturally  in  terms  of  S  and  P  waves;  each  mode  is  repre¬ 
sented  by  a  Fourier  coefficient.  This  allows  the  energy  contribution  from  each  mode 
to  be  calculated.  These  quantities  will  be  useful  to  determine  the  magnitude  of  P  to  S 
conversion,  and  vice  versa,  and  which  features  of  the  medium  dictate  the  magnitude 
of  the  conversion.  Included  in  this  discussion  is  the  matrix  formulation  of  the  method, 
and  a  detailed  description  of  the  numerical  procedure  for  making  the  recursion  matrix 
unitary. 


In  chapter  3,  we  test  the  applicability  of  the  method  for  a  2D  sample  problem 
whose  exact  solution  is  also  calculated.  This  provides  a  check  on  the  method  and  some 
guidelines  to  follow  concerning  the  range  of  parameters  for  which  the  approximation 
is  valid.  This  will  help  in  the  extrapolation  of  the  phase-screen  method  to  problems 
whose  exact  solutions  cannot  be  obtained.  The  actual  comparison  is  made  in  two  ways. 
First,  by  plotting  the  cartesian  components  of  both  solutions  for  various  parameters 
of  the  problem,  and  second,  by  computing  the  average  of  the  norm  of  the  difference 
of  the  two  wave  vectors.  To  perform  these  comparisons  more  efficiently,  we  first 
determine  the  marginal  number  of  phase  screens  needed  such  that  adding  additional 
phase  screens  for  fixed  observation  distance  does  not  significantly  improve  the  answer. 

Chapter  3  is  concluded  with  an  examination  of  energy  flux  conversion  for 
the  test  problem.  Using  a  single  phase  screen,  we  extremize  the  S  to  P  and  P  to 
S  conversion  as  a  function  of  the  ratio  of  the  width  of  the  regions  to  the  free  space 

2  In  deed  the  limit  as  the  distance  between  phase  screens  approaches  *ero  has  some  resemblence  to 
the  path  integral  approach  to  wave  propagation,  neglecting  backscattering,  Dashen  [1979]  and  Flattl 
[1983]. 


wavelength  for  various  ratios  of  the  two  sound  speeds.  This  wiih  indicate  winch  wave¬ 
lengths,  or  equivalently  what  structure  sizes,  induce  the  most  conversion,  and  how  it 
depends  on  the  velocities. 


In  chapter  4,  the  phase-screen  method  for  random  media  is  discussed,  and 
the  details  of  -some  specific  examples  are  given.  We  are  primarily  interested  in  the  rate 
of  conversion  when  the  medium  is  characterized  by  random  velocity  fluctuations,  and 
how  it  depends  on  the  frequency,  the  correlation  length,  ^e  velocity  fluctuations,  the 
distance  of  propagation,  and  the  type  of  initial  disturbance.  In  particular,  we  model  a 
two  dimensional  random  medium  with  only  one  length  scale  by  randomly  aligning  the 
phase  screens  and  ensemble  averaging  over  realizations.  Each  realization  is  specified 
by  a  set  of  computer  generated  psuedo-random  numbers  that  determine  by  how  much 
each  phase  screen  gets  shifted.  We  find  that  the  rate  of  energy  conversion  has  a  robust 
dependence  on  the  magnitude  of  the  velocity  fluctuations  and  the  structure  size,  and 
to  a  lesser  extent  on  the  frequency  and  whether  the  initial  disturbance  is  a  P  wave  or 
an  S  wave.  We  also  find  that  the  energies  tend  to  constant  equilibrium  values  after 
propagating  a  distance  far  enough,  depending  on  the  parameters  which  determine 
the  rate  of  conversion,  and  that  the  equilibrium  ratio  of  S  wave  to  P  wave  energies 
is  roughly  given  by  the  ratio  of  the  P  wave  to  S  wave  phase  velocities.  Although 
this  is  merely  a  precursor  of  the  actual  3D  problem  with  multiple  length  scales  to  be 
investigated,  it  demonstrates  the  applicability  of  the  method  to  address  questions  of 
discrimination. 


And  finally,  in  chapter  5  we  draw  some  conclusions  from  our  exercises,  and 
discuss  work  currently  in  progress,  which  will  utilize  the  the  techniques  developed 
in  this  report  to  address  more  pertinent  and  realistic  problems  related  to  seismic 
discrimination. 
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SECTION  2 


THE  METHOD 


2.1  DERIVATION  OF  THE  RECURSION  RELATIONS 


The  phase-screen  method  may  be  formulated  as  an  initial  value  problem. 
Between  any  consecutive  pair  of  phase  screens  the  uniform  elastodynamic  wave  equa¬ 
tion  for  the  displacement  vector,  using  the  average  values  for  the  parameters  of  the 
medium,  is  solved.  The  unique  solution  in  that  region  is  determined  by  specifying 
the  initial  value  of  the  displacement  on  the  first  of  the  two  screen.  This  initial  data 
is  obtained  by  evaluating  the  displacement  in  the  previous  region  at  this  screen  and 
modifying  the  phase  by  a  position  dependent  phase  factor.  The  solution  in  each  region 
may  be  expressed  as  a  Fourier  expansion,  decomposed  into  S  and  P  waves.  The  Fourier 
coefficients  of  the  successive  solution  may  be  determined  in  terms  of  the  Fourier  co¬ 
efficients  of  the  previous  solution  and  the  phase  factor.  Thus  by  knowing  the  initial 
displacement  produced  by  a  source,  the  displacement  after  N  phase  screens  may  be 
recursively  determined. 

The  elastic  wave  equation  for  the  displacement  vector  in  a  uniform  medium, 
and  in  the  absence  of  external  forces  is 

=  fiV’u  +  (A  +  M)V{V  ■  u),  (10) 

where  u  is  the  displacement  vector,  /i  and  A  are  the  average  Lam4  moduli,  and  p  is  the 
average  density  of  the  medium.  The  linearity  of  this  equation  allows  the  solutions  to 
be  decomposed  into  S  and  P  waves.  Thus  the  displacement  vector  may  be  expressed 
as  the  sum 


u  =  Up  +  Us, 


(11) 


where  the  two  pieces  are  chosen  to  satisfy  the  constraints 

Vxu^O  (12) 
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V  •  us  =  0. 


(13) 


It  may  readily  be  shown  that  this  separation  can  always  be  accomplished,  and  that  it 
is  unique.  Using  this  decomposition  the  wave  equation  reduces  to  two  simple  homo¬ 
geneous  wave  equations  of  the  form 


1  d*u  P 
c*P  at* 


=  VV 


(14) 


1  d*us 

4  at * 


=  V*us, 


(IS) 


where  ej»  and  cs  are  the  mean  propagation  speeds  for  the  respective  disturbances, 
given  in  terms  of  the  density  and  Lam6  moduli  by 


(16) 

(17) 


The  general  form  of  the  solution  of  the  these  equations  that  propagates  in 
the  forward  direction  is 


* 


u(r,l)  =  J  <tt,it,[«„4(t.,*,)ei‘'-'  +  l|fl.(t.,i,).‘k>’je-i*',  (18) 

where  A{k9tk¥)  and  Ba(kw ,kp)  are  Fourier  coefficients  for  the  P  wave  and  the  two 
polarisations  of  the  S  wave  respectively,  and  r  =  (*,*).  There  is  an  implicit  sum 
over  a  from  1  to  2  in  this  expression,  following  the  standard  convention  of  summing 
over  repeated  indices  unless  otherwise  specified.  For  simplicity,  the  factor  of  e'**4 
will  be  suppressed  throughout  the  rest  of  the  analysis  since  a!)  of  the  solutions  we  are 
interested  in  here  will  have  the  same  time  dependence.  Note  that  the  three  dimensional 
Fourier  expansion,  which  is  typically  used  to  express  solutions  of  this  type,  hu  been 
reduced  to  a  double  integral  by  use  of  the  dispersion  relations,  which  allow  the  wave 
vectors  to  be  expressed  solely  In  terms  of  ka  and  h,  as 
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kp =(*„*„(< »’/4  - 


(19) 


is  =  (*..*,.  (<»74  -  *i  -  (20) 

The  longitudinal  and  transverse  unit  vectors,  Ip  and  l§,  are  defined  such  that  the 
constraints  of  eqs.  (12,13)  are  satisfied,  ».  e.  kp  x  Ip  =  0,  and  k$  •  eg  =  0,  a  =  1,2. 
The  unit  vectors  satisfying  these  constraints  are 


ep  =  kp  =  —{ktfk„kpt) 


i's  =  £*4  =  -o  -  ^)-v*(t,„o,-t.) 
\jxks\  w1 


(its  x  IJI  w*  to1  Cg  v 


where  kp,  and  Jfcs,  are  the  ^-components  of  eqs.  (19,20)  respectively. 

For  a  given  set  of  Fourier  coefficients,  propagation  between  the  phase  screens 
is  accomplished  using  eq.  (18).  Let  us  now  introduce  a  phase  screen  at  z  =  A*  to 
compensate  for  the  previously  ignored  detailed  structure.  Schematically,  if  the  initial 
displacement  vector,  u^,  is  specified,  then  the  resultant  displacement  vector  beyond 
the  screen,  denoted  by  u^,  is 


(2<i 

The  superscripts  on  the  u’s  will  be  used  to  indicate  the  number  of  phase 
screens  the  wave  has  been  propagated  through.  The  solution  for  any  value  of  the 
superscript  takes  on  a  slightly  modified  form  of  the  general  solution  of  eq.  (18).  After 
N  phase  screens  the  displacement  vector  is 

nW(r)  .  I fkfaAWtknkJ +»*-*** 

+  IS  (*«.**)  €*••**»*.  (28) 
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The  Fourier  coefficients  have  been  defined  with  the  exponential  factors  that  depend 
on  N  already  factored  out  to  simplify  the  initial  value  equation  at  z  =  NAz. 

At  this  point  we  depart  from  the  traditional  phase-screen  approximation  for 
scalar  waves,  which  has  been  realised,  as  in  eq.  (24),  by  multiplying  the  entire  wave  of 
the  previous  region  by  an  overall  position  dependent  (but  independent  of  ka  and  ifc,) 
phase  to  provide  the  initial  data  for  the  solution  in  the  next  region.  In  our  approach 
we  will  multiply  the  longitudinal  and  transverse  contributions  by  different  A-phases, 
which  depend  on  the  appropriate  phase  velocity  for  each. 

Furthermore,  each  Fourier  mode  is  multiplied  by  its  own  phase  factor, 
i.  e.  the  A-phases  depend  on  ka  and  ka.  This  simply  reflects  the  fact  that  the  phase 
the  wave  accumulates  from  propagating  some  distance  in  the  z-direction  depends  on 
the  angle,  relative  to  the  z-axis,  at  which  it  propagates.  This  angle  is  given  by 

arctan((*J  +  Aj)l/*/(^*/c*  ~*l~  *!)1/1)*  (26) 

For  small  angle  propagation  the  phase  is,  however,  independent  of  the  angle  to  first 
order,  and  hence  independent  of  ka  and  ky.  This  is  satisfied  if  kT,  kv  <  w/c ,  such  that 
k,  cs  w/c.  This  will  be  true  for  either  scalar  or  vector  waves  if  the  typical  length  scales, 
in  the  plane  perpindicular  to  the  z*axis,  are  much  greater  than  w/c  ( c  represents  either 
the  transverse  or  longitudinal  speed  for  the  case  of  vector  waves).  We  are,  however, 
interested  in  problems  of  elastic  wave  propagation  with  length  scales,  frequencies,  and 
speeds  that  do  not  strictly  satisfy  this  condition,  and  hence  the  phase  accumulated  is 
not  independent  of  ka  and  fcy. 

Incorporating  these  two  extra  features,  if  A^~l^[kail t¥)  and  Bls~^(katkv) 
are  the  coefficients  of  the  expansion  of  ».  e.  the  wave  having  been  already 

propagated  through  JV  -  1  phase  screens,  then  the  initial  value  equation  for  may 

be  expressed  as 


+  «}  j  ,***“»',  (27) 

where  the  uniform  phase  that  was  acquired  in  propagating  from  (N  -  l)Az  to  N Az 
has  been  absorbed  into  the  definitions  of  the  A-phases. 
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The  coefficients  A^(k)  and  flJW(k)  may  be  determined  in  terms  of  A^_1)(k) 
and  B[N~l^(k)  and  the  A-phases  by  setting  the  RHS  of  this  equation  equal  to  the  RHS 
of  eq.  (25),  evaluated  at  z  —  NAz,  and  then  multiplying  both  sides  of  this  new  equa¬ 
tion  by  e-'***“'k»v,  integrating  over  x  and  y  from  -oo  to  +oo,  and  projecting  onto  the 
appropriate  polarization  vectors  defined  above.  The  result  of  this  procedure  is 


ks'kp  J  * 


=  r— — . —  [itdyl d'k'  X  i°s(k'„k',))  • 

ks'kpj  t 

■  |(fc,  x 

+  (i,  x  (») 

where  the  vector  products  may  be  determined  from  cqs.  (19-23). 

As  a  quick  check  that  the  results  of  eqs,  (28,29)  are  sensil*?^  ssAs  that  if  th^ 
A-phases  are  independent  of  z  and  y,  then  the  z  -  y  integration  may  be  performed 
to  yield  a  product  of  delta  functions  of  V,  -  k,  and  k*  ~  k*.  The  integral  over  d1# 
may  also  be  performed  now.  replacing  k'  everywhere  with  fc.  Evaluating  the  vector 
products,  it  becomes  dear  that  the  contribution  to  involves  only  and  not 

The  analogous  result  is  true  for  each  of  the  as  well.  This  shows  that  if 
the  medium  is  truly  uniform,  there  is  no  interc  on  version  between  the  P  and  two  types 
of  S  waves.  However,  if  the  A-phases  do  depend  on  the  *  -  y  coordinates,  as  they 
will  for  interesting  problems,  the  coefficients  of  the  P  wave  after  the  N**  phase  screen 
depends  on  both  polarisations  of  the  S  wave  and  the  P  wave  of  the  preceding  region; 
a  similar  result  is  true  for  the  S  waves  u  well.  Thus  the  A-phases  encode  information 
about  the  scattering  and  wave  type  conversion  as  the  wave  propagates  through  an 
inhomogeneous  medium. 


13 


The  A-phnses  at  the  N °*  phase  screen  in  terms  of  the  velocity  anomalies 
between  {N  —  1)A z  and  NAt  are 


»  > 


«/j 


(30) 


and  a  similar  expression  with  P— »  S.  At  this  point  the  phase-screen  method  has  been 
reduced  to  finding  expressions  for  6cp(x,  y,z')  and  £es(x,y,z')  for  the  physical  prob¬ 
lem,  and  specifying  the  initial  form  of  the  wave.  In  the  next  chapter  we  will  look  at 
a  particular  example  where  their  spatial  distribution  is  given  deterministically  from 
the  medium.  Eventually,  however,  they  will  be  treated  as  stochastic  quantities.  First, 
however,  there  are  some  further  general  results  to  be  obtained. 


3.2  MARGINAL  NUMBER  OF  PHASE  SCREENS 


To  determine  the  number  of  phase  screens  to  be  used,  a  self-contained 
method  is  needed.  In  the  next  chapter  we  will  have  an  exact  solution  with  which 
to  compare.  We  could  increase  the  number  of  phase  screens  per  observation  distance 
until  the  two  solutions  agree  as  well  as  they  will.  For  most  problems,  however,  the 
exact  solution  is  not  accessible.  We  will  determine  the  marginal  number  of  phase 
screens  needed  (MNPS  for  short)  such  that  the  phase-screen  solution  converges  to 
its  asymptotic  answer  to  within  some  given  tolerance  as  the  number  of  phase  screens 
are  increased.  This  is  self-contained  and  will  optimise  computer  time,  using  only  the 
necessary  number  of  phase  screens  that  yields  an  answer  close  to  the  the  asymptotic 
answer  one  would  obtain  if  there  were  an  infinite  number  of  phase  screens. 

More  explicitly,  to  determine  the  MNPS,  a  Cauchy  criterion  of  sorts  is  used. 
The  distributions  of  velocity  anomalies  are  chosen,  and  then  the  phase-screen  answer 
is  computed,  adding  more  and  more  phase  screens  for  a  fixed  observation  distance. 
The  MNPS  is  determined  when  the  norm  squared  of  the  difference  between  succes¬ 
sive  displacement  vectors,  one  computed  with  K  more  phase  screens  than  the  other, 
averaged  over  x  and  y,  converges  to  within  a  specified  tolerance. 

For  convenience,  we  write  the  displacement  vector  as 

U'"l(r)  =  /  {»«'«>(*..*,..)  +  )«!“>(*,,  *,,i)  + 

(3*1 
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where  the  superscript  M  in  this  expression  denotes  the  total  number  of  phase  screens 
used  for  a  given  observation  distance.  The  norm  squared  of  the  difference  of  two 
displacement  vectors,  both  the  same  except  t*>at  one  is  calculated  using  K  more  phase 
screens  than  the  other,  is 


=  II  ji'k{i  («!“>(*)  -«!"-*>(*)) 

+  }  („<«><*)  -  .««-*>(*))  +  i  («<*>(*)  -  |p 

=  f  t i‘k  d'e  { (»!">(*)  - 

+  (“{“’(‘I  -  uj “-*>(*))  (o '“>(*')•  -  „<“-*)(*')■) 

+  («<«>{*)  -  »<«-*>(*))  (»'«>(*')•  -  „!«-*!(*•)•)  j. ‘I*. 

(32) 


Integrating  this  quantity  over  x  and  y,  and  dividing  by  the  length  of  the  interval  in 
each  direction,  yields  a  product  of  delta  functions  of  the  wave  numbers  as  the  lengths 
of  the  intervals  go  to  infinity.  The  average  of  the  norm  of  the  squared  difference, 
denoted  ANSD,  is  given  by 


ANSD{s/)  ^ 

+  l<4V)(*»*/)  ~  */)!*} 

<  tolerance.  (33) 


The  smallest  value  of  M  such  that  this  criterion  is  satisfied  will  determine  the  number 
of  phase  scree  uj  to  be  used. 
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2.3  ENERGY  FLUX 


lb  conclude  the  general  results,  we  compute  an  expression  for  the  time* 
averaged  energy  flux  for  both  the  longitudinal  and  transverse  displacements  in  terms 
of  the  Fourier  coefficients  of  the  phase-screen  method.  The  energy  flux  will  also  be 
averaged  over  the  x  -  y  plane,  which  allows  the  total  energy  to  be  expressed  as  the 
sum  of  the  individual  energies  of  the  P  and  the  two  components  of  S  waves.  Thus 
as  a  function  of  s,  the  ratios  of  P  to  S  and  S#  to  Sy  (the  two  polarisation  of  the  S 
wave)  may  be  simply  determined.  Also,  the  ratio  of  transmitted  to  incident  fluxes 
may  be  extremised  with  respect  to  the  parameters  of  the  medium,  providing  a  means 
to  determine  what  types  of  regions  produce  significant  or  small  wave  type  conversion. 

The  energy  flux  given  in  terms  of  the  displacement  is 


-A(Vu) 


dui 

li 


(34) 


where  u  in  this  expression  represents  only  the  real  part  of  the  complex  displacement, 
and  the  subscripts  refer  to  the  cartesian  components.  Time-averaging  this  expression 
produces 


Si  «!*{-*(*..) ~  » E  +  |")  (,un*P J  •  (W) 

where  now  u  represents  the  full  complex  displacement  without  the  time-dependent 
exponential  phase. 

Inserting  the  displacement,  given  by  eq.  (25),  into  this  expression,  and  av¬ 
eraging  over  the  x  -  y  plane,  the  only  nonvanishing  component  of  the  averaged  fluxes 


are  the  r-components,  which  for  the  Af* 

region  may  be  expressed  as 

(30) 

82 Ls*  1 

(35) 
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where  a  -  1,2  and  \k\  -  +  Jk*.  The  limits  of  integration  are  due  to  the  fact 

that  for  larger  )£)  the  inodes  are  evanescent,  and  do  not  contribute  a  real  part  to  the 
flux,  as  should  be  expected  for  a  non-dissipative  medium.  From  these  expression!!  it  is 
obvious  that  energy  is  conserved  throughout  the  region  between  any  consecutive  pair 
of  phase  screens.  Energy  conservation  from  one  region  to  the  next  is  not  nearly  so 
trivial,  and  wilt  now  be  examined. 


2.4  MATRIX  FORMULATION  AND  ENERGY  CONSERVATION 


The  examination  of  energy  conservation  is  most  efficiently  performed  using 
a  matrix  representation  of  the  recursion  relations.  Although  there  exists  a  analogous 
matrix  representation  in  three  dimensions,  and  a  generalization  for  continuous  values 
of  the  wave  number,  for  the  purposes  of  the  rest  of  this  report,  and  to  clarify  the  details, 
we  shall  consider  two  dimensional  media  with  periodic  boundary  conditions  in  the  x- 
direction.  The  2D  problem  may  be  recovered  as  a  special  case  of  the  3D  problem  by 
setting  kv  «  0,  and  the  periodic  boundary  conditions  replace  the  continuous  variable 
kt  with  mv/2d,  where  the  period  is  id  and  m  is  an  integer;  the  corresponding  integrals 
are  replaced  with  sums.  The  Fourier  coefficients  BXm  decouple  from  A.  and  BXm  hi 
the  recursion  relations  when  kv  «  0,  and  hence  will  be  set  to  zero. 

Dropping  the  ‘1*  subscript  from  B\m,  we  shall  rescale  the  Fourier  coefficients 
as 


m 


m 


and  define 


where  Afj» 
tivdy,  and 


the  vector  of  Fourier  coefficients 


(Ai*i  .(*s  Am  Aw  Am  M{s)  «(.v) 


.5 


w 

I 


(«) 

aad  M,  we  (he  greatest  iatagen  teaa  than  Tiv/ncr  and  gdw/ae,  mtpec- 
T  denotes  the  transpose  of  the  row  vector.  Only  the  modes  with  harmonk 


IT 


e-dependence  are  included  in  this  vector.  For  integers  larger  than  these  values  the 
modes  are  evanescent,  and  will  not  be  indu  vcd.  This  vector  may  be  written  schemat¬ 
ically  in  block  form  as 


(^\ 

IV’1"1)  =  .  («) 

\  Bi">  I 

and  the  matrix  that  represents  the  recursion  relations  may  also  be  expressed  in  block 
form  as 

I C  v”\ 

<«) 

V"  V’J) 

where  the  integers  m,n  range  from  -Afj»  to  +Af^,  and  the  integers  ji,i/  range  from 
-Ms  to  +M5.  The  olf-diagonal  blocks  represent  conversion  of  one  type  of  wave  into 
another;  their  explicit  form  will  be  determined  shortly.  The  recursion  relations  for  the 
Fourier  coefficients  may  now  be  expressed  as 

(  \  (U”  ViZ 

UH  Wi'  vs 

where  the  summation  convention  is  employed,  or  equivalently  as 


r  l  > 


\V^  a>U\Vt"m%  (44) 

In  addition  to  setting  i,s0  and  making  i*  discrete,  the  A-phascs  may  be 
expanded  in  Fourier  series  as 

B  ^  pM  §***•&,  (45) 

* 

and  a  similar  expression  for  the  S  wave  4-phase.  Combining  this  information  with 
eqs.  (25*29)  yields 
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The  total  energy  Uux  in  the  region,  aver  aged  over  time  and  the  x~ 
coordinate,  is 

$\S)  »  W{FW|V(A0).  (50) 

2 

This  simple  form  wa*  made  possible  by  the  rescalings  of  the  Fourier  coefficients  in 
eqs.  (38,39).  The  energy  of  the  P  wave  is  given  by  the  first  2 My  +  1  terms  of  this 
expression,  and  the  energy  of  the  S  wave  is  given  by  the  remaining  terms.  Using  the 
recursion  relations  this  may  be  re-expressed  as 

a  (51) 

Thus  fur  energy  to  be  conserved  from  one  region  to  the  next  the  recursion  matrix  U 
must  be  unitary,  t.  e.  if  /  is  defined  to  be  the  identity  matrix,  then 

V'V  »  UV'  «  /.  (S3) 
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For  general  A-phases  the  recursion  matrix  is  not  unitary,  although  for  the 
range  of  parameters  for  which  the  phase-screen  method  is  expected  to  be  valid,  it  is 
almost  unitary.  Almost  unitary  is  defined  to  mean  that  the  eigenvalues  of  the  matrix 
are  close  to  being  unimodular,  and  the  commutator  of  the  matrix  with  its  adjoint  has 
small  entries.  Although  the  violation  of  energy  conservation  at  each  phase  screen  may 
be  less  than  a  tenth  of  a  percent,  the  cumulative  energy  violation  after  propagating 
through  a  thousand  phase  screens  may  be  significant.  To  remedy  this  problem  we 
have  employed  the  following  algorithm. 

First,  the  complex  vector  |F)  is  written  as  a  real  vector  of  twice  its  original 
length  simply  by  entering  the  real  and  imaginary  part  of  each  complex  entry  as  pairs  of 
entries  in  the  real  vector.  The  complex  matrix  U  is  also  put  in  real  form  by  replacing 
each  complex  entry  Uij  with  the  2  x  2  block 


l  ReUij  —ImUij  ' 
K  ImUij  ReUij  ; 


(S3) 


If  the  real  matrix  corresponding  to  U  is  denoted  by  R,  the  condition  for  energy  to  be 
conserved  is  that  R  be  an  orthogonal  matrix,  defined  to  satisfy 


RrR  =  RRt  =  /. 


(54) 


The  second  step  of  the  algorithm  is  to  perform  a  singular  value  decomposi¬ 
tion  of  Ry  i.  t.  write 


R  =  OxDOt  (55) 

where  Ox  and  02  are  orthogonal  matrices,  and  D  is  a  diagonal  matrix  whose  entries  are 
the  singular  values  of  R.  The  singular  values  of  a  matrix  are  the  positive  square  roots 
of  the  eigenvalues  of  RrR.  This  decomposition  must  be  performed  by  a  numerical 
routine  in  practice.  The  matrix  R  can  be  said  to  be  almost  orthogonal  if  its  singular 
values  are  all  close  to  unity.  To  make  R  exactly  orthogonal  we  simply  replace  D  with 
the  identity  matrix  /.  Thus 


R  — ►  0i02, 


(56) 
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which  is  clearly  orthogonal  since  Oj  and  O2  are  both  orthogonal. 

At  this  point  we  have  the  necessary  machinery  to  work  on  specific  problems, 
which  is  what  we  will  turn  to  now. 


SECTION* 


COMPARISON  WITH  AN  EXACT  SOLUTION 


We  shall  illustrate  the  phase-screen  method  for  a  problem  whose  exact  so¬ 
lution  is  tractable  without  being  trivial.  The  test  problem  we  have  chosen,  shown  in 
figure  3,  is  to  solve  for  the  elastic  displacement  vector  in  a  two  dimensional  medium 
constructed  from  two  different  homogeneous  strips  of  equal  width,  repeating  alter¬ 
nately  and  infinitely  in  the  ^-direction,  and  of  infinite  extent  in  the  2-direction,  taken 
to  be  the  forward  direction  of  wave  propagation.  We  shall  center  the  origin  at  the 
mid-point  of  region  one,  and  by  convention  let  the  regions  both  be  a  distance  2d  wide. 


The  exact  solution  for  this  problem  is  calculated  first,  and  then  the  phase- 
screen  approximation  is  computed  and  compared  to  it.  We  shall  use  this  comparison 
to  quantitatively  understand  the  limitations  of  the  method.  In  particular,  we  examine 
the  limits  of  the  region  widths,  the  magnitudes  of  velocity  fluctuations  about  the 
mean,  and  the  propagation  distances  for  which  the  approximation  is  valid.  We  also 
provide  some  explanation  for  why  these  limits  exist.  This  analysis  will  yield  some 
guidelines  for  the  range  of  parameters  over  which  the  phase-screen  method  may  be 
trusted  for  problems  whose  exact  solutions  are  not  accessible. 

Also,  the  ratio  of  the  flux  for  a  transmitted  S  wave  to  the  flux  of  an  incident 
P  wave  is  extremized,  for  a  single  phase  screen,  with  respect  to  the  ratio  of  the  width 
of  the  regions  to  the  free  space  wavelength.  This  analysis  is  performed  for  various 
ratios  of  cp/cs-  The  same  result  is  computed  for  an  incident  S  wave  and  transmitted 
P  wave.  This  result  yields  insight  into  the  wavelengths,  or  equivalently  structure  sizes, 
that  contribute  the  most  to  wave  type  conversion,  and  the  dependence  on  the  wave 
speeds. 


3.1  THE  EXACT  SOLUTION 

We  begin  by  presenting  the  exact  solution  to  the  problem  stated  above. 
By  “exact”  we  mean  that  there  is  no  fundamental  approximation  made  for  the  wave 
equation  or  the  boundary  conditions  used  to  solve  this  problem.  Some  numerical  work 
must  be  performed,  however,  to  obtain  the  final  result,  which  cannot  be  written  m 
dosed  form,  but  the  solution  may  be  determined  to  any  desired  accuracy. 
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The  {procedure  (for  solving  this  problem  is  relatively  straightforward.  .First, 
the  homogeneous  .wave  equation  is  solved  in  each  region.  Then  these  solutions  are 
matched  together,  determining  the  unique  rotation,  by  satisfying  the  boundary  con* 
ditionsat  tthe  interfaces  and  theinitial.conditionfor  the  wave  at  *  =  0. 

To  simplify  finding  the  solution,  the  symmetries  of  the  problem  may  be  ex¬ 
ploited.  First,  themediwn  possesses  periodic  translation  invariance  in  the  x-direction. 
Any  ;pair  Of  neighboring  regions,  labelled  tone  and  two,  are  indistinguishable  {torn  all 
other  such  pairs.  Only  the  solution  tin  a  (fundamental  .pair  of  regions  must  be  calcu¬ 
lated;  .the  ‘boundary  conditions  are  matched  only  at  fhe  interlace  between  these  two 
regions.  Thcadlutiona  in  all  other  pairs  of  regions  are  identical  to  this  one. 

Second,  the  medium  is  dearly  invariant  under  reflections  about  the  2-axis. 
Solutions  may  (be  separated  into  parity  eigenfunctions  that  are  either  even  or  odd  as 
x  -*  —x.  If  even  and  odd  parity  solutions  exist  and  may  be  treated  separately  in  each 
region,  the  number  of  unknowns  to  tolve  for  in  the  boundary  value  equations  will  be 
reduced  by  a  factor  of  two.  Clearly  the  solution  in  the  region  centered  on  the  2-axis, 
from  x  =  -d  to  x  =  d,  may  be  separated  into  even  and  odd  parity  eigenfunctions.  We 
shall  .now  argue  that  the  solution  In  all  of  the  other  regions  must  be  of  the  same  parity, 
when  expanded  about  the  midpoint  Of  that  particular  region,  as  the  one  centered  on 
the  x-axis. 

Refer ing  to  figure  3,  parity  about  the  z-axis  alone  only  guarantees,  for  ex¬ 
ample,  that  the  solution  between  x  =  d  and  x  -  3d  is  an  even  or  odd  reflection  of 
the  solution  between  x  --  -3d  and  x  =s  -d.  However,  if  we  simultaneously  invoke 
translation  invariance,  there  is  also  parity  invariance  about  ail  of  the  lines  parallel  to 
-the  2-axis  and  centered  at  the  midpoint  of  each  region.  Furthermore,  the  parity  of  the 
solution  about  these  lines  must  be  the  same,  either  even  or  odd,  in  all  of  the  regions, 
for  the  solutions  to  match  smoothly  at  the  interfaces. 

To  see  that  these  statements  are  true  we  argue  that  the  convene  results 
leads  'to  contradictions.  Suppose  the  solutions  in  the  regions  from  x  »  d  to  x  ==  3d 
and  from  x  -  -3d  to  x  @  -  d  respect  the  parity  symmetry  about  the  s-axie,  but 
not  about  the  tines  centered  at  the  midpoints  of  the  two  regions.  It  is  trivial  to  see 
that ’these  solutions  cannot  be  identical.  This  contradicts  the  asumption  that  there  is 
translation  invariance.  Pushing  this  argument  slightly  further,  suppose  the  solution 
in  the  region  from  x  -  -d  to  x  -  d  has  a  particular  parity  about  the  s-axis,  and 
the  solution  in  the  region  from  x  -  d  to  x  -  3d  has  the  opposite  parity  about  its 
.midpoint.  The  solution  in  the  region  from  x  «  -3d  to  *  s  ~d  that  is  identical  to  the 
notation  in  the  region  from  x«dtox»3d  cannot  join  smoothly  to  the  solution  in 
region  one  at  x  =  -d,  i.  c.  the  derivative  of  the  displacement  with  respect  to  x  will 
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be  discontinuous  there.  Thus  the  boundary  condition  involving  such  derivatives  will 
not  be  satisfied,  and  hence  this  is  not  a  valid  solution.  These  arguments  are  trivially 
extended  to  the  other  regions  as  well 

Combining  these  symmetries,  we  may  arbitrarily  choose  the  fundamental 
pair  of  regions  to  be  the  one  between  x  =  -d  and  x  =  3d;  the  solutions  in  these 
two  regions  may  be  separated,  both  into  either  even  or  odd  eigenfunctions,  expanded 
about  the  midpoint  of  that  region.  For  example,  the  x-dcpendcnce  of  the  solution 
between  x-d  and  z  =  3d  may  be  expressed  in  terms  of  either  even  or  odd  functions 
of  the  argument  (z  -  2d).  The  boundary  conditions  at  the  interface  of  these  two 
regions  have  thus  been  greatly  simplified. 

The  wave  equation  to  solve  is  the  same  as  eq.  (10)  except  that  here  there  are 
two  regions,  each  with  their  own  densities,  clastic  moduli,  and  hence  wave  velocities. 
In  each  region  the  general  solution  corresponds  to  that  of  a  homogeneous  displacement 
wave.  The  general  solution  in  region  one  may  be  expressed  as 


U.M  =  Jdk,  '  '] ,  (57) 

and  to  obtain  the  solution  in  region  two  simply  let  1  2  and  x  (z  -  2d). 

To  solve  for  the  coefficients  of  the  expansion  it  is  necessary  to  match  the 
boundary  conditions  at  the  interface.  It  is  worth  noting  that  we  have  written  the 
expression  above  in  terms  of  k,  instead  of  k,  as  for  the  phase-screen  case  to  simplify 
matching  the  boundary  conditions  for  all  values  of  z.  Using  the  dispersion  relations, 
the  wave  vectors  in  each  region  may  be  written  in  terms  of  k,  as 


(58) 

k*  =  ((*74  -  k\)'i\k.). 

(59) 

As  discussed,  the  solutions  may  be  separated  by  their  parity.  Note  that  in 
order  to  satisfy  the  constraints  of  eqs.  (13, 13)  if  the  x-component  of  the  displacement 
vector  has  a  particular  parity,  the  component  must  have  the  opposite.  Thus  there  are 
only  two  parity  eigenfunctions  that  the  solution  may  decomposed  into,  which  for  the 
sake  of  convenience  will  be  labelled  (even, odd)  or  (odd, even)  according  to  whether  the 
z-component  has  even  or  odd  parity  respectively.  Note  also  that  the  initial  condition 
chosen  must  be  consistent  with  the  parity  of  the  solution,  c.  y.  if  the  incoming  wave 
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is  a  uniform  plane  P  wave,  the  solution  should  have  even  parity  for  the  component, 
and  if  the  incoming  wave  h  a  uniform  plane  8  wave,  the  solution  should  have  even 
parity  for  the  x-component.  We  will  examine  the  latter  case  first.  For  simplicity,  the 
following  definitions  will  be  used: 


jl/2 

(60) 

»,-at 

tt> 

(61) 

cp,k 
a  =  — 
w 

(62) 

(  cUJ\ 

(63) 

—H 

% 

where  the  subscript  refers  to  region  one  or  two.  Using  these  definitions  the  solution 
in  region  one  with  (even, odd)  parity  may  be  written  as 


ut(r)  53  ^dfcjjiaicosfwajx/c*.,)  +  icji'sinfwajx/cj*,)]  A\{k) 

+  [i6|CO«(wdtx/c5,)  -  idjisinfwdjx/e^)]  fij(4)|ea*.  (64) 

There  is  a  similar  expression  in  region  two  with  1  -*  2  and  x  (x  -  2d).  All  of  the 
square  roots  denote  the  principal  value  of  the  complex  square  root,  and  k  >  0  since 
only  the  forward  wave  is  included. 

Now  we  will  proceed  to  determine  the  unknown  coefficients  by  solving  the 
boundary  conditions  at  x  =  d.  There  are  two  vector  boundary  conditions  at  the 
interface  to  be  satisfied.  First,  the  displacement  must  be  continuous  at  the  interface 
to  insure  that  the  velocities  are  finite  everywhere.  And  second,  the  traction  (the 
normal  projections  of  the  stress  tensor)  must  be  continuous  there  ss  well  to  avoid 
infinit#  acceleration*. 

Continuity  of  displacement,  U|(d_)  =  us(d*),  yields  the  following  two  equa¬ 
tions: 
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aj  cos  ( waid/cpt )  A\(k)  +  61  cos  (tudid/cs,)  B\(k)  — 

aj  cos  ( wa^d/cp, }  Ai(k)  +  6j  cos  (wdjd/cs^)  Bj(fc)  (65) 


ci  sin  ( wa\djcpx )  Ai(k)  -  d\  sin  (wdidfcs,)  B\(k)  — 

-  cjsin  (u/ajd/cp,)  At(k)  +  dj  sin  (todjd/csJ  Bi(k).  (66) 

Similarly,  continuity  of  traction  at  *  =  d,  Tj‘(d_)  =  2?(d+),  yields  two  equations. 
First,  the  normal  projection  of  the  stress  tensor  may  be  written  in  terms  of  the  dis¬ 
placement  as 

t-'.-Wv.  «-,.(§£ +  §■£).  (67) 

To  clarify  the  notation,  all  superscripts  here  refer  to  cartesian  components  of  tensors 
or  vectors,  and  xl  =  x,  x1  =  a.  The  subscripts,  as  for  the  displacement  vector,  will 
refer  to  the  particular  region.  Using  the  expression  for  the  displacement,  eq.  (64),  the 
two  traction  equations  are 


Hx  I  ™ axCxCos(waxdjCf^)Ai(k )  +  — (b\  -  d|)  SQs{v>dxd/eSl)Bx(k)  > 

l  <?.  **  1  . 

Pi  I  —aiCjcos(u?oid/cp,)Ai(fc)  +  ~t~(6j  -  «4)cos(u/d2d/cs,)I?i(k)  |  (68) 

l  Cs*  * 


-i-  fzpio?  +  Aj|  5in(u>aid/ep,)4i(l:)  +  Pidthi8in(u>did/es,)fl|(k)  — 

ept  1  *  cst 

-  «  [2pjo|  +  Aj|sm(waid/ep,)4*(t)  -  ^^jds6ism(u>djd/es,)0i(*)*  (69) 

There  are  four  sets  of  unknown  Fourier  coefficients,  and  one  set  of  discrete 
wave  numbers  to  solve  for.  We  will  divide  the  four  boundary  condition  equations 
arbitrarily  by  flj(k),  use  three  of  the  equations  to  solve  for  Bj B»,  and  At/ B%, 
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and  then  insert  these  expressions  into  the  fourth  equation  to  determine  the  discrete 
eigenvalues  allowed  for  kt  labelled  km,  where  m  is  an  integer.  The  remaining  set  of 
discrete  coefficients,  {BJm}  =  {Bt(km)},  specifies  the  shape  and  overall  normalisation 
of  the  displacement  wave,  which  is  given  by  an  initial  condition  that  we  will  choose  at 
a  =  0. 


As  will  be  shown  shortly,  the  eigenvalue  equation  is  a  transcendental  equa¬ 
tion,  making  it  necessary  to  solve  for  the  eigenvalues  numerically.  This  forces  us  to 
choose  specific  values  for  the  parameters  of  the  problem.  A  choke  that  simplifies  the 
algebra  considerably  is  to  let  Aj  =  Aj  =  0  and  /ij  =  /ij.  The  densities  of  the  two 
regions  must  be  different  to  have  an  interesting  problem,  but  by  making  this  choke 
for  the  elastic  moduli,  their  explicit  dependence  cancels  out  of  the  transcendental 
equations,  and  only  the  phase  velocities  need  to  be  given  numerical  values. 

The  four  equations  are  divided  by  Bj,  however,  the  solutions  of  the  first 
three  equations  for  Aj/Bj,  Bt/Bj,  and  Aj/Bj  are  ail  divided  by  the  determinant  of 
the  matrix  of  Cramer’s  rule  when  solving  a  linear  system  of  equations.  We  will  absorb 
this  determinant  into  the  definition  of  Bj(Jb),  and  drop  the  subscript  so  that  we  may 
write  the  solution  in  a  symmetric  form  for  the  two  regions.  Using  the  first  three 
boundary  conditions,  along  with  the  choices  for  the  Lame  moduli  just  mentioned, 
some  tedious  but  straightforward  algebra  yields  the  following  expressions: 


=  Cj(6j  -  &j)sin(woad/c#«,)cas(ti/djd/cs,)coa(ti/djd/c$,) 


4-  aj6jdj sin(wdjd/c^) cos(wajd/c^)cos(wdid/ej,) 


4-  08^jdisin(u>d|d/c(s,)co8(wo!d/c^)cos(wdjd/es,)  (70) 


T\(k)  m  B\(k)fB{k)  s  «J^el(un(wo|d/c^)cos(wasd/c4i)cos(wdid/e5l) 

4-  o2kte)$in(tt>08d/c^)cos(wa)(l/c^i)cos(wd|d/es*)i  (71) 

and  similar  expressions  for  £|  and  7t,  but  with  1  *-*  2.  7*  is  simply  the  determinant 
of  the  matrix  discussed  above. 

Inserting  these  expressions  into  the  fourth  equation,  the  transcendental 
eigenvalue  equation  for  k  is 


0  a  tti\(waid/cpi)t^ii(wdid/esi)axbiCidi  +  i*n{waid/ef>i)tm(wdid/est)aibiC\di 
+  tan(u>oj  d/ept)  Unlwdid/cs^a^cjdi  +  tan(u>aid/ej>,)  tan(u>djd/ cs,)aj6jcjdj 
+  tan(u;a1d/cpl)Un(tyojd/e^)(6j|-4*)(cI-ei).  (72) 

It  has  been  written  in  terms  of  the  tangent  function  by  dividing  by  the  necessary 
sines  and  cosines.  Writing  it  in  this  form  simplifies  the  numerical  root  finding  pro¬ 
cess  because  the  asymptotes  of  the  tangent  functions  provide  natural  limits  to  search 
between  for  the  roots.  In  fact,  between  any  consecutive  pair  of  asymptotes  there  are 
either  two  roots,  one  root,  or  no  roots. 

All  that  remains  to  completely  specify  the  solution  is  to  fix  the  initial  shape 
and  normalization  of  the  wave,  specified  by  the  set  of  coefficients  {Bm}.  To  do  this 
we  will  choose  that  at  z  =  0 


ut(s:,0)  =  u,(i,0)  =  *.  (73) 

This  expression  yields  four  equations,  but  there  is  only  one  set  of  coefficients  {£«}.  In 
general,  each  of  these  four  conditions  must  be  satisfied  by  a  unique  linear  combination 
of  complete  functions.  Each  linear  combination  b  specified  by  the  set  of  coefficients 
that  gets  multiplied  to  these  complete  functions.  Thb  implies  that  unless  the  func¬ 
tions,  that  have  been  used  to  expand  the  solution,  are  equivalent  to  four  complete 
sets  of  functions,  there  b  no  hope  of  satisfying  the  four  initial  condition  equations  at 
s  5=  0;  the  solution  will  be  over-determined.  It  b  worth  noting  that  the  cosine  and 
sine  functions  used  to  expand  the  solution  are  not  orthogonal  functions,  making  it 
impossible  to  determine  by  casual  inspection  the  number  of  complete  set#  of  functions 
they  comprise.  However,  an  examination  of  the  eigenvalues  of  k  reveab  that  they  may 
be  grouped  into  pairs.  Thb  b  much  more  apparent  for  the  larger  values  of  it.  In  fact, 
the  discrete  pairs  asymptotically  approach  either  inx/d  or  »(n  +  t/2)x/d.  Combining 
thb  with  the  fact  that  there  are  both  sine  and  cosine  functions  present  in  the  solution 
indicates  that  there  are  actually  enough  complete  sets  of  functions  to  satisfy  the  initial 
condition.  Although  thb  b  a  only  a  heuristic  argument,  we  can  proceed  to  actually 
solve  the  initial  condition  to  verify  that  a  solution  does  exbt. 

Another  consequence  of  the  fact  that  the  cosine  and  sine  functions,  used 
to  expand  the  solution,  are  not  orthogonal,  b  that  projecting  the  solution  for  the 
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displacement  onto  the  initial  condition,  a*  is  done  in  typical  Fourier  analysis,  does  not 
yield  the  solution  for  the  individual  coefficients.  Instead  we  solve  for  the 
numerically  using  a  least  squares  routine.  The  variable  x  b  first  dberetbed  to  the  set 
{*»),  and  the  sum  over  the  discrete  values  of  fc^is  truncated.  Also,  the  set  {£*}  may 
be  written  as  a  column  vector  of  dimension  M,  where  M  is  the  value  of  m  at  which 
the  series  b  truncated.  The  initial  condition  may  also  be  written  as  a  column  vector 
of  dimension  N%  where  N  b  four  times  the  number  of  discrete  values  of  x  chosen, 
since  there  are  two  spatial  vector  components  in  each  of  the  two  regions.  The  linear 
transformation  that  maps  the  coefficients  into  the  initial  condition  b  an  N  x  M  matrix 
whose  entries  depend  on  km  and  x».  To  solve  for  the  coefficients  requires  inverting 
thb  matrix,  which  b  performed  numerically  by  the  least  squares  routine.  To  make 
thb  dearer,  schematically  the  initial  condition  at  z  -  0  may  be  written  as 

*l\  l  ' 

v.sums  =  °  .  AO 

\*uJm  \0  )m 

where  the  subscript  n  labels  the  discrete  values  of  x.  Thus  each  entry  actually  repre¬ 
sents  a  column  of  N/4  entries.  Also,  the  general  solution  at  t  »  0,  may  be  written  in 
the  form 

(75) 

m 

The  solution  for  the  unknown  coefficients  {$„}  b 

B- =  £«««.•  (76) 

m 

Having  solved  this  inversion  numerically,  and  using  eqs.  (64,70,71),  the  com¬ 
plete  solution  for  the  displacement  for  -d  <  x  <  d  with  (even, odd)  parity  b 


*i(r)  =  £|[fo|*€os(w«»»x/«#i)  +  iei.ssui(a>«i.x/e^)]£i» 

+  Jfo»»eos(wd|»x/ej,)  -  £d,..i*in(wd*.*/e*,)]  Tj.Jfl.e'4*.  (77) 


and  the  solution  for  d  <  x  <  3d  Is  the  simitar  except  with  1  —  2  and  z  (z  ~  2d). 

Having  calculated  the  solution  above,  it  is  now  straightforward  to  write 
down  the  (odd, even)  parity  solution.  In  region  one,  the  displacement  is 

Ui(r)  “  / dfc | (ioji sin  (w«|X/ej>J  +  kcx  cos  (u>a,x/c*>j]  Ax(k) 

+  |thjt  sin  (wdxz/cs,)  -  kdx  cos  (wd»x/ej,)]  Bx (A)Je<k*,  (78) 

and  as  before,  a  similar  expression  for  region  two  with  1  -*  2  and  x  (z  -  2d). 

Continuity  of  displacement  at  x  —  d,  Uj(d_)  «  Uj(d+),  yields  the  following 
two  equations: 

aj  sin  (watd/cpt)  A,  (k)  +  &jsin (wdxd/cst)  Bx(k)  = 

-  a:  sin  (v.  z9df  cft)  A,(k)  -  5,  sin  (wd3d/ e*J  Jtt(k)  (79) 


cj  cos  (waxd/cpt)  4|(fc)  -  dj  cos  (wdxdjegt)  Bx[k)  a 

cj  cos  (u/ajd/cnj)  4* (A)  -  dt  cos  (u>d:d/ e$J  B3(k) .  (80) 

Similarly,  continuity  of  traction  at  x  «  d,  TJl{d_)  ®  Tjl(d,),  yields 


Mt  {  -^e|Ctsln(watd/ej04i(A)  +  -  d|)#in(ttHM/c*J0i(*})  = 

-  Mt  {^o,e,sin(w«|d/e^)4|(b)  +  ~(h|  -  df )  sia(tt>d,d/cs,)lM*)  J  (81) 


~  (2pi«!  +  A»]  cos(we»d/e*)4t(*)  +  ^pid|8|Cos(wd>d/c*,)ll|(A)  « 

—  fepjof  +  A,]  ce*(watd/e*)4|(fc)  4  -3'iiidt**ct*(wdtd/cjJ#t(A).  (82) 
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Am  before,  using  the  fint  three  boundary  condition  equation*,  choosing  At  = 
At  a  0  and  mi  =  Mr,  end  absorbing  the  determinant  of  Cramer’s  rule  into  the  definition 
of  B(A),  the  coefficients  are 

Ci[k)  =  Ai[k)/B[k)  -  ct(4|  -  ^)cos(w«ad/c^)sin(i»d|d/cs,)ain(«»dj</c*l) 

+  «sltdtCos(u>djd/csa)sin(iMsd/ci^)ain(«Mf|d/cjt) 

+  *aMi  co^a>dtd/c5k)  sin(wa,d/e*,)  sinfwdtd/c*)  (83) 

T»(fc)  5  Bi(k)/B{k)  =  a,&,ejCO«(u><i,d/c^)  sinfuKiid/c*,)  sinfwdjd/c*,) 

+  Mlbtctcm(waid/t^)»iu{wald/e^)sin(w4td/cst)t  (84) 

and  similar  expressions  for  £s  and  T*,  but  with  1  2. 

Inserting  these  expressions  into  the  fourth  equation,  the  transcendental 
eigenvalue  equation  for  A  is 


0  s  tan(w«|d/c|«()Un(wdtd/cjv)a|^e|di  +  tan(wo)d/c^)tan(wd|d/es))«|A|e)dt 
+  tan(w«sd/c^)  tan(wdtrf/cs,)ethiejdt  -¥  taa(wo:<i/c^)  taa(w4|d/e*ja*AiCtdt 
♦  tan(wd|d/csj  tan(w4^/cij)(*5  -  Aj)(e{  - 1\).  (85) 

The  initial  condition  at  «  s  0  tor  the  (odd, even)  parity  case  is 

Ut(x,0)  B  U|(z,0)  a  A,  (86) 

and  as  before,  the  coefficients  will  be  solved  for  using  a  least  squares  routine.  With 
these  coefficients,  the  solution  may  be  written  as 
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Ui(r)  “  E{ [«W«n  (wai»*/ei»()  +  kcXm cos (walmzf cPl ) j  £lm 

+  (i&i*isin  (wdx *i/cj,)  -idim cos (wdi,*x/«s,)]  (87) 

For  any  other  initial  condition  that  is  uniform  as  a  function  of  x,  the  solution 
will  be  a  linear  combination  of  the  displacement  vectors  given  by  eqs.  (77,87).  For 
initial  conditions  that  do  depend  on  x,  the  only  change  in  the  solution  is  that  the  set 
of  coefficients  {Bm)  must  be  determined  in  the  same  manner  as  before,  but  for  the 
new  initial  condition.  This  completes  the  analysis  of  the  exact  solution,  and  we  now 
turn  to  the  phase-screen  solution  of  the  same  problem. 


3.2  THE  PHASE  SCREEN  SOLUTION 


Now  that  the  exact  solution  has  been  calculated,  we  would  like  to  compute 
the  phase-screen  approximation  for  the  same  physical  problem,  and  then  compare  the 
two.  Figure  4  depicts  how  the  problem  will  be  treated  using  the  phase-screen  method. 
The  wave  is  propagated  uniformly  between  phase  screens,  as  discussed  in  chapter  2. 
The  phase  of  the  wave  at  the  phase  screen  is  then  retarded  uniformly  for  values  of 
x  corresponding  to  region  one,  and  advanced  uniformly  for  values  of  x  corresponding 
to  region  two.  Because  of  the  periodicity  of  the  problem,  the  general  solution  of  the 
uniform  wave  equation,  written  as  a  Fourier  integral  in  chapter  2,  may  be  expressed 
as  a  discrete  series.  Thus  the  displacement  vector  between  the  Nih  and  the  (JV  +  l)** 
phase  screens  may  be  written  as 

um(r)  =■  + 

m 

+  (hU  -  &w)  }  e<w/w,  (88) 

where  are  of  the  same  form  as  those  defined  in  eqs.  (60,61,62,63)  except 

that  the  k  refers  to  the  x-component  of  the  wave  vector  in  this  case,  given  by  k 
m*/2d,  and  instead  of  carrying  a  suhsript  that  depends  on  the  wave  velocities  of  the 
two  regions,  there  is  only  one  quantity  of  each  type  that  depends  on  the  average  of 
either  the  Iran  verse  or  longitudinal  velocities  of  the  two  regions. 

At  this  point  the  velocity  anomalio  need  to  be  specified  for  the  problem. 
For  the  I*  wave  in  region  one 
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-3d  -d  +d  +3d 


{♦—  region  1  -*[*-  region  2  “*] 


Figure  4.  The  phase-screen  treatment  of  the  test  problem  introduces 
the  phase  screens,  represented  by  the  heavy  solid  lines.  The 
square  wave  shows  schematically  how  the  phase  is  retarded  or 
advanced  at  the  phase  screens. 
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cp  +  6cp(x,  z)  =  epu  (89) 

and  in  region  two  let  1  -*■  2  in  this  expression.  Using  eq.  (30)  the  argument  of  the 
A-phase  for  the  P  wave  is 


(®j)  J  A z  if  -d,  <  x  <  d 

)  J  Az  if  d  <  x  <  Zd 


(90) 


This  is  actually  a  periodic  function,  which  simply  repeats  for  the  other  regions.  We 
shall  express  the  A-phase  as  the  Fourier  series 


n 


where  for  n  =  0 


n(p)  _  I 
Vm0  ~  2 


1/2 


A  z 


(92) 


and  for  n  /  0 


£><£>  =  - 
“  2 


-<»)*) 


1/2 


Am 


sin  nn/2 

m r/2 


(93) 


There  are  also  similar  results  for  the  S  wave  with  P  S  in  the  expressions  above. 

Inserting  the  A-phases  into  the  expressions  derived  in  eqs.  (46-49)  produces 
the  recursion  matrix  for  this  problem.  Recall  that  the  Fourier  coefficients  in  the  matrix 
formulation  have  been  rescaled  according  to  (38)  and  (39) .  In  terms  of  the  notation 
used  here,  the  elements  of  this  matrix  are 


4?  =  A.,o,  (99) 

and  for  the  (odd, even)  parity  initial  condition  of  eq.  (86),  let  *-*  in  these  two 
expressions. 

As  discussed  for  the  exact  solution,  for  any  other  initial  conditions  inde¬ 
pendent  of  x,  linear  superpositions  of  these  expressions  may  be  used  because  of  the 
linearity  of  the  wave  equation.  For  initial  conditions  that  do  depend  on  z,  elementary 
Fourier  analysis  may  be  used  to  determine  the  initial  coefficients. 

In  principle,  the  phase  screen  solution  for  the  medium  we  are  interested  in 
is  complete.  To  actually  obtain  the  displacement  after  N  phase  screens  requires  nu¬ 
merical  work  to  iterate  the  recursion  relations  N  times.  We  will  numerically  compare 
the  phase-screen  solution  to  the  exact  solution  in  the  next  section. 
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3.3  COMPARISON  OF  THE  SOLUTIONS 


To  compare  the  two  solutions  requires  numerical  routines  to  obtain  both 
solutions,  and  a  method  of  comparison.  Although  there  are  many  ways  we  could 
compare  the  two  answers,  we  have  chosen  to  do  so  in  the  following  two  ways.  First, 
the  cartesian  components  of  the  two  displacement  vectors  at  fixed  values  of  x  are 
plotted.  Also,  the  average  difference  of  the  vectors  at  fixed  z  is  computed.  Explicitly, 
this  means  that  the  absolute  norm  of  the  difference  of  the  vectors  is  computed  at  a 
given  value  of  z,  integrated  over  z,  and  divided  by  the  length  of  the  interval  to  obtain 
the  average  difference. 

Before  proceeding  with  these  comparisons,  the  marginal  number  of  phase 
screens,  discussed  in  chapter  two,  must  be  calculated.  For  this  periodic  problem  the 
Cauchy  criterion  derived  in  eq.  (33)  becomes 

ANSD(Zf)  -  [d  z/)  - 

=  £  {kmH*/)  -  +  !*£?(*/)  “  «im"X)(*/)l*} 

m 

<  tolerance ,  (100) 

where 

t»M(«/)  =  cmAWe<wa-,'/<,’(M+,)  +  </w^)e<wd-^/es(w+1)  (101) 

u^}(*/)  =  (102) 

Note  that  for  M  phase  screens  evenly  distributed  over  a  distance  z/,  Az  ~  z//(M  + 1). 
We  have  chosen  the  tolerance  to  be  0.005  for  the  following  results,  and  incremented 
the  number  of  phase  screens  by  adding  10  for  every  50  km  of  observation  distance, 
t.  e.  K  =  0.2z/  is  used  in  eq.  (100),  Tables  1  and  2  contain  the  MNPS  results  for 
a  series  of  observation  distances  using  various  combinations  of  the  frequency,  width, 
and  magnitude  of  velocity  fluctuations. 

From  these  tables  a  general  guideline  for  the  number  of  phase  screens  to  be 
used  may  be  extracted; 
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Table  3. 


Marginal  number  of  phase  screens  for  c$i  -  3.6  km/s, 
est  =  3.8  km/s,  cp/cs  =  y/2  . 


tv  (rad/s) 

*/  (km) 
d  (km) 

50 

100 

400 

5.0 

10.0 

100 

200 

400 

2.5 

10.0 

160 

320 

5.0 

200 

400 

5.0 

320 

800 

10.0 

2.5 

800 

1600 

Marginal  number  of  phase  screens  for  tSs  -  3.5  km/s, 
e$,  s=  3.9  km/s,  tries  -  >/2  • 


Hi 

■1 

too 

0 

400 

wmmm 

10.0 

170 

340 

1360 

10.0 

5.0 

340 

680 

1360 

3360 

2.S 

10.0 

240 

480 

960 

2560 

5.0 

5.0 

480 

960 

2560 

5120 

10.0 

2.5 

320 

640 

1600 

3200 
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where  L  is  used  here  to  represent  the  observation  distance,  a  —  2  x  10s  is  the  constant 
of  proportionality  for  the  tolerance  we  have  chosen,  and  0  ~  2.  For  a  less  stringent 
tolerance  the  number  of  phase  screens  needed  will  be  much  less. 

Using  this  information,  the  comparison  is  made  as  efficiently  as  possible.  The 
parameter  space  to  be  examined  includes  the  magnitude  of  the  velocity  fluctuations 
and  three  length  scales,  the  wavelength,  the  width  of  the  regions,  and  the  observation 
distance.  One  of  these  length  scales  may  be  eliminated,  however,  by  noting  that  the 
solution  is  invariant  under  rescaling  all  of  the  lengths  by  the  same  factor.  Thus  we 
are  free  to  fix  the  wavelength  and  vary  the  remaining  parameters. 

Figures  5-16  contain  plots  of  the  real  parts  of  the  cartesian  components 
of  the  phase-screen  and  exact  displacement  vectors,  both  in  the  same  figures,  at  a 
sequence  of  observation  distances.  The  solid  lines  correspond  to  the  exact  solution 
while  the  dashed  lines  correspond  to  the  phase-screen  solution;  the  x-component  of 
the  displacement  is  plotted  above  the  r-component.  Also,  the  average  of  the  difference 
of  the  norm  of  the  two  vector  solutions  has  been  included  in  the  figure  captions.  For 
all  of  these  comparisons  ep/cs  =  v/2,  which  was  chosen  to  simplify  the  analysis  of  the 
exact  solution. 

These  comparisons  demonstrate  the  accuracy  and  qualitative  features  of  the 
method,  although  we  have  included  only  a  small  sampling.  For  example,  it  is  clear 
that  the  accuracy  of  the  method  decreases  for  larger  observation  distances  and  smaller 
widths,  as  compared  with  the  wavelength,  and  for  larger  velocity  perturbations.  Fur¬ 
ther  analysis  shows  that  the  quantitative  conditions  for  the  method  to  be  valid  are 

kd>  I,  (104) 

such  that  backscattenng  is  small,  and 


where  again  0  -  2  for  the  parameters  we  are  using.  This  latter  condition  is  very 
similar  to  a  geometrical  optics  condition,  and  is  demonstrated  by  comparing  figures  S- 
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Figure  6.  Comparison  of  the  exact  and  pbatMcran  solutions  for 

*>  =  »d/*i  c*  a  3.T  km/«,  4  =  S.0  km,  =  0.1  km/s, 

*/  a  100.0  km,  eve.  diff.  =  0.041. 


*/<3 

Figaro  T.  Compvtaoo  of  tho  «uct  and  phju*-*cr«*n  KdutioM  for 

ir  =  5.0  rad/t,  I*  =  5.T  km/*,  4  =  5.0  km.  =  0.1  km/r, 
«/  a  200X)  km,  avo.  diff.  s  0.052. 


u((a 


*  =  5.0  r*d/m  g,  =  5.7  km/*,  4  a  54)  km,  U$  =  0.1  km/*, 
*/  »  400  0  km,  ave.  diif.  =  0.141* 
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ux(x)  vs.  x/d 


u*(x)  vs  */d 
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Wp*»  Comparboo  of  tha  mad  and  phm  icrwa  mIuIIou  for 

*  =  *  0  rad/*,  it  =  3.7  km/»,  d  =  10.0  km,  U$  =  0.1  km/*, 
*/  =  500  km,  ava.  dift  =  0410. 
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•/a 

Figure  10.  Comparison  of  lb*  exact  and  phase*#r«een  solutions  for 

w  s  5.0  ra d/s,  if  s  5.T  km/«,  4  s  10.0  km*  4c#  =  0.1  km/s. 
a/  »  100.0  km,  ave.  diff.  =  0.024. 


i 
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•/* 

Flgurt  11.  Comparison  of  tba  exact  and  phatMcrem  solutions  for 

m  s  i.O  rad/s,  f«  =  J.T  km/s,  4  =  10.0  km.  Icj  =  0.1  km/s, 
«/  *  200.0  km.  avt.  dlff.  =  0.030. 


Figure  13.  Camper  Uoo  of  the  met  end  ph**fr-*ere«m  solutions  for 

tv  --  5.0  r»d/s*  i$  =  3.7  ktu/s,  4  ss  10.0  km*  U$  =  0.1  km/s, 
«/  *  4000  km,  eve.  diff.  =  0.040. 


ux(x)  VS.  x/d 


uz(x)  vs.  x/d 


Figure  13.  Comparison  of  the  exact  and  phase*  screen  solutions  for 

w  =  5.0  rad/s,  =  3.T  km/s,  d  —  10.0  km,  =  0.2  km/s, 
zj  -  50.0  km,  ave.  diff.  =  0.035. 
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x/d 

Figure  14.  Comparison  of  the  exact  and  phase-screen  solutions  for 

w  —  5.0  rad/s,  c$  =  3.7  km/s,  d  =  10.0  km,  6cs  =  0.2  km/s 
zj  -  100.0  km,  ave.  diff.  =  0.050. 


Figure  15.  Comparison  of  the  exact  and  phase-screen  solutions  for 

w  as  54)  rad/s,  eg  —  3.7  km/a,  d  =  10.0  km,  6c s  =  0.3  km/s, 
Zj  w  200D  km,  ave.  diflf.  =  04)66. 


iparison  of  the  exact  and  phase-screen  solutions  for 
5.0  rad/s,  cs  =  3.7  k m/s,  d  =  10.0  km,  6cs  —  0.3  km/s, 
400.0  km,  ave.  dlff.  =  0.114. 
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condition,  and  is  demonstrated  by  comparing  figures  5-8  to  9-12.  Since  the  width  in 
the  latter  figures  was  doubled,  the  displacement  may  be  propagated  to  four  times  the 
distance  at  the  same  accuracy.  Also,  comparing  figures  9-12  to  13-16  demonstrates 
that  that  doubling  the  velocity  fiuctuation  reduces  the  observation  distance  by  a  factor 
of  four  to  maintain  the  same  accuracy.  These  relations  have  been  verified  by  other 
similar  comparisons. 


3.4  ENERGY  CONVERSION 


We  conclude  this  chapter  with  an  examination  of  the  parameters  that  induce 
the  maximum  conversion  of  energy  from  one  type  of  wave  to  another.  To  accomplish 
this  we  maximize  the  energy  conversion  across  a  single  phase  screen  for  an  incident 
plane  wave.  The  dependence  of  conversion  on  the  velocity  fluctuations  comes  from  the 
Fourier  coefficients  of  the  A-phases,  given  in  eqs.  (92,93);  the  absolute  value  squared  of 
these  coefficients  enter  into  the  expression  for  the  energy.  These  factors  are  bounded 
from  above  regardless  of  the  magnitude  of  the  velocity  perturbations,  and  satisfy 


|A*>1*  <  1  (106) 

IA-.I’  <  (107) 

wr 

The  dependence  on  the  velocity  perturbations  is  eliminated  by  using  the  combination 
of  the  perturbations  and  A z  that  give  the  maximum  conversion.  Using  the  recursion 
matrix,  given  by  eqs.  (94-97),  with  the  optimum  phase  amplitude  on  the  screen,  the 
conversion  is  maximized  for  various  ratios  of  cp/cs  with  respect  to  the  dimensionless 
variable 


X  = 


2dw 
i rep  ’ 


(108) 


where  the  factor  of  2  and  x  are  due  to  the  convention  of  making  each  region  a  distance 
2d  wide. 


Table  3  contains  the  values  of  x*  labelled  Xm**>  such  that  the  conversion 
is  maximized  for  various  ratios  of  the  sound  speeds.  It  illustrates  that  the  P  wave 
converts  energy  at  a  faster  rate  than  the  S  wave  except  when  the  sound  speeds  are 
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equal,  and  that  the  length  scales  that  induce  the  most  conversion  of  either  type  are 
approximately  equal  to  the  wavelength  of  the  P  wave.  This  regime  is  considered 
“resonance  scattering”,  and  contributes  the  most  to  scattering  effects;  for  x  <  1  the 
medium  is  quasi-homogeneous,  white  for  x  >  1  the  medium  becomes  transparent, 
Wu  (1988).  The  entries  of  Xma*  =  1-00  in  table  3(a)  are  due  to  the  fact  that  the 
conversion  is  maximized  as  this  dimensionless  variable  gets  smaller,  but  for  x  <  1<0, 
Mp  =  0  (cf.  section  2.4),  and  the  recursion  matrix  becomes  trivial,  not  admitting 
conversion.  Finally,  note  that  as  the  ratio  cp/cs  increases,  the  P  wave  loses  energy 
more  rapidly,  while  the  S  wave  loses  energy  more  slowly. 


Table  3. 

(a) 


Maximum  energy  conversion  from  one  phase  screen  for  an  initial 
(a)  P-wave,  (b)  S-wave. 


(b) 


Cp/Cs 

1.0 

V2 

2.0 

4.0 

Xmas 

1.17 

1.00 

1.00 

1.00 

sfi/sg 

0.31 

0.44 

0.58 

0.64 

cp/cs 

1.0 

v/2 

2.0 

4.0 

Xmas 

1.17 

1.11 

1.06 

1.01 

4!i /sJ8 

0.31 

0.24 

0.21 

0.19 
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SECTION  4 

RANDOM  MEDIA 


Here  we  apply  the  phase  ecr«n  method  to  media  for  which  the  position  de¬ 
pendent  velocity  anomalies,  te #>(x,  y,  a)  and  fc,j(x,  y,  a),  are  characterised  statistic  ally. 
They  will  be  treated  as  sero  mean  random  variables.  Eventually  we  will  employ  phase 
screens  with  multiple  wavelengths  to  describe  random  inhomogeneities  on  a  variety  of 
length  scales.  The  random  phases  for  this  problem  may  be  computed  in  terms  of  the 
power  spectral  density  (PSD)  of  the  medium,  Knepp  {1963]  and  Martin  and  Flatt£ 
(1988),  but  this  will  not  be  done  here.  Instead  we  will  limit  ourselves  to  one  length 
scale  to  obtain  the  results  that  follow. 

To  model  a  random  medium  with  only  one  characteristic  length  scale  we 
use  the  same  periodic  phase  screens  as  were  used  in  the  previous  chapter,  aligned  over 
some  correlation  length,  but  then  randomly  shifted  by  some  distance  less  than  the 
wavelength  of  the  screen  in  the  x-direction,  over  the  next  correlation  length.  This 
process  is  then  repeated  to  propagate  the  wave  to  the  observation  point  of  interest. 
The  number  of  phase  screens  used  per  propagation  distance  for  the  random  medium 
problem  will  be  the  same  as  the  MNPS  determined  for  the  layered  problem,  for  par¬ 
ticular  values  of  the  widths  of  the  regions,  frequencies,  mean  velocities,  etc. 

In  general,  a  random  shift  in  the  alignment  of  the  phase  screen  relative  to 
the  positioning  of  the  original  one  is  accomplished  by  letting 


A»(x)  — ♦  Aw(x  -  4 <h),  (109) 

where  is  a  random  number  such  that  0  <  *y  <  1.  Refering  back  to  eqs.  (91-97), 
the  affect  of  this  shift  on  the  recursion  matrix  is  straightforward  to  determine;  the 
elements  of  the  recursion  matrix  transform  according  to 


(110) 


(HI) 


(112) 


USJ  — >  U” 


USJ  —  t£s  (113) 

where  again  ~Mp  <  m,n  <  Mp  and  -Ms  <  n,u  <  Ms- 

The  transformation  that  implements  the  random  shift  may  be  represented 
as  a  matrix  as  well,  written  in  block  form  as 


(  0  \ 

r  =  .  (114) 

V  0  e1*^  J 

Note  that  only  the  diagonal  entries  arc  nonvanishing,  and  there  is  no  implicit  sum¬ 
mation  in  this  expression.  The  transformation  law  for  the  recursion  matrix  under  a 
random  shift  is  then  given  by 


v  — >  ri/r 1  (us) 

It  is  easily  verified  that  T  is  a  unitary  matrix,  or  orthogonal  when  expressed 
in  its  equivalent  real  form.  Thus  we  need  only  modify  the  recursion  matrix  U  once 
to  make  it  unitary;  all  random  transformations  of  the  modified  recursion  matrix  will 
automatically  lead  to  total  energy  conservation. 

For  the  simulations  we  perform,  the  recursion  matrix  of  the  problem  in 
chapter  3,  given  by  eqs.  (92-97),  will  be  employed.  The  energy  flux,  averaged  over  x 
and  i,  is  computed  as  a  function  of  observation  distance  for  a  number  of  realisations. 
Each  realisation  is  specified  by  a  set  of  computer  generated  psuedo- random  numbers 
which  determine  the  alignment  of  the  set  of  phase  screens.  We  shall  then  ensemble 
average  over  realizations  to  smooth  out  spurious  results.  In  practice,  we  shall  average 
over  one  hundred  such  realisations. 

Figures  17-21  plot  the  averaged  energy  flux  as  a  function  of  distance  out  to 
2000  km,  and  for  a  variety  of  parameters.  Although  we  are  not  attempting  to  model 
the  lithosphere  here,  the  length  scales  and  velocity  perturbations  are  roughly  those 
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found  in  the  lithosphere  at  NORSAR,  Wu  and  Aki  {1085],  where  the  correlation  length 
is  10  -  20  km  and  the  rms  velocity  perturbations  are  2%~4%.  The  total  energy  flux 
has  been  normalised  to  one  by  dividing  by  the  incident  flux.  The  pairs  of  figures  on 
each  page  are  produced  for  identical  media,  including  the  particular  realizations  that 
are  used;  the  upper  figures  are  for  initial  P  waves,  while  the  lower  figures  are  for  initial 

I  S  waves. 

i 

The  significant  features  of  energy  conversion  in  this  2D  random  medium  are 
that  the  conversion  rate 

r 

1.  for  P  to  S  occurs  at  a  faster  rate  than  for  S  to  P; 

2.  is  approximately  independent  of  the  frequency  for  high  frequencies; 

3.  is  inversely  proportional  to  the  length  scale  of  the  medium;  and 

4.  is  proportional  to  the  square  of  the  velocity  fluctuations. 


These  results  may  be  qualitatively  compared  to  the  Born  approximation 
results  obtained  by  Wu  and  Aki  (1985].  The  first  result  is  simply  due  to  the  fact  that 
the  impedance  of  the  P  wave  is  greater  than  that  of  the  S  wave,  and  the  last  three  are 
typical  of  high  frequency  scattering.  In  particular,  the  lack  of  frequency  dependence 
is  easily  seen  by  noting  that  the  displacement  scales  inversely  with  the  frequency,  but 
the  energy  flux  is  given  by  the  frequency  squared  times  the  displacement  squared, 
which  does  not  scale  with  the  frequency. 

Kindly,  the  energy  plots,  particularly  figures  20  and  21,  illustrate  that  the 
energies  of  the  waves  tend  toward  fixed  values  if  propagated  far  enough.  These  values 
are  independent  of  the  initial  conditions,  and  the  asymptotic  ratio  of  the  S  wave  to  P 
wave  energy  is  roughly  e #»/««. 

For  propagation  in  a  3D  random  medium,  realised  in  the  way  described 
above,  the  expected  result  is  that  the  energy  will  be  equally  partitioned  among  the  S 
wave  polarisations  and  weighted  according  to  the  ratio  of  the  phase  velocities  relative 
to  the  P  wave. 
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f LUX  fRANS/tNC  FLUX  TRANS/inC 


INITIAL  P-WAVE 


INITIAL  S-WAVE 


Figure  17.  Energy  conversion  for  a  random  medium  with  w  s  $.0  rad/s, 
d  —  S.O  km,  i$  =  3.7  km/*,  Se$  =  0.1  km/*,  e*/e*  =»  \/i. 


5 


Figure  tl  &wcy  couwiIod  for  *  random  modlmn  with  w  —  1041  rtd/i, 
4  =  8.0  km,  I|  s  IT  km/t,  8c«  s  04  km/c,  c^/tf  *  4. 
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Figure  19.  Energy  conversion  for  •  random  medium  with  w  =  S.0  n 
d  =  10.0  km,  is  =  3*T  km/t,  =  0.1  km/t,  c*fcs  «  >/i 


MMN5/INC  FLUX  TRAWS/lNC 


INITIAL  P-WAVE 


INITIAL  S-WAVE 


Figort  90.  Eurgy  wwiiot  for  «  random  medium  with  «  s  LQ  ntd/s. 
d  SS  9j0  km.  Ms  ss  9.7  km/6*  fe#  =  0.3  km/*. 
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31.  Energy  e<mv«r«ba  fbr  a  random  medium  with  w  =  *.o  rad/*, 

4  s  5,0  km,  i$  m  M  km/*,  #<$  as  0.3  km/e,  r?U9  ®  3* 

01 

1 

...  I 


SECTION  S 


CONCLUSIONS  AND  FUTURE  WORK 


In  this  report  we  have  developed  the  phase-screen  method  for  elastic  vec¬ 
tor  waves,  and  demonstrated  its  use  to  approximate  elastic  wave  propagation  in  two 
dimensions.  By  comparing  with  an  exact  solution  we  were  able  to  determine  the  gen¬ 
eral  conditions  for  which  the  method  is  valid.  The  problem  in  chapter  3  proved  to 
be  quite  demanding,  particularly  on  the  number  of  phase  screens  needed  to  approxi¬ 
mate  the  exact  solution.  This  was  due  to  the  fact  that  the  layered  structure  allowed 
only  a  small  number  of  the  modes  to  propagate  for  parameters  of  interest,  inducing 
substantial  scattering  at  oblique  angles  relative  to  the  direction  of  propagation  of  the 
original  wave.  For  smoother  phase  distributions,  e.  g.  Gaussian  phase  screens  etc.,  we 
expect  that  the  method  will  require  fewer  phase  screens,  and  may  reduce  to  the  usual 
“parabolic  approximation”  to  produce  accurate  results. 

We  also  demonstrated  how  the  method  may  be  used  to  compute  P  to  S 
conversion  in  terms  of  the  averaged  energy  flux.  We  found  that  it  is  energetically 
more  favorable  for  a  P  wave  to  convert  its  energy  into  an  S  wave  at  a  faster  rate  than 
for  the  reverse  process,  particularly  when  the  longitudinal  phase  velocity  is  much 
larger  than  the  transverse  phase  velocity.  Furthermore,  it  was  demonstrated  that  the 
greatest  conversion  occurs  for  structure  sizes  comparable  to  the  P  wave  wavelength. 

The  application  of  the  method  to  a  random  medium  was  initiated  by  ran¬ 
domly  aligning  the  phase  screens  and  ensemble-averaging  the  results.  This  provides 
a  genuine  first  step  towards  treating  a  truly  random  medium  problem.  The  results  of 
the  analysis  here  showed  how  the  energy  conversion  again  proceeds  at  a  faster  rate  for 
an  initial  P  wave  until  the  energies  reach  a  fixed  asymptotic  value.  The  rate  of  conver¬ 
sion  is  proportional  to  the  magnitude  of  the  velocity  fluctuations  squared,  inversely 
proportional  to  the  structure  size,  and  roughly  independent  of  the  frequency  for  high 
frequencies.  Also,  the  asymptotic  ratio  of  S  wave  to  P  wave  energies  is  approximately 
given  by  c p/c$.  These  results  indicate  that  using  P-S  ratios  as  a  discriminant  is  vi¬ 
able  on  regional  distance  scales  of  up  to  2000  km,  provided  the  velocity  perturbations 
are  not  too  large. 

Although  this  work  demonstrates  many  of  the  strengths  of  the  phase  screen 
method,  it  also  indicates  the  need  for  a  great  deal  of  further  work  to  be  able  to 
tackle  more  realistic  problems.  The  following  features  to  be  incorporated  into  the 
phase-screen  framework  are  under  current  investigation: 


First,  the  3D  random  medium  problem  with  one  characteristic  length  scale  is 
being  developed.  This  will  allow  S#  —  Sy  ratios  as  well  as  P— S  ratios  to  be  computed. 
Second,  the  formulation  in  the  time-domain  is  being  developed.  This  will  allow  the 
behavior  of  transient  solutions  to  be  investigated,  and  the  time  delay  separating  the 
P  wave-train  from  the  S  wave-train  to  be  computed.  Third,  we  are  interested  in  using 
sources  with  spherical  symmetry  to  model  explosions.  And  fourth,  the  incorporation 
of  phase  screens  with  all  length  scales  represented,  weighted  according  to  the  PSD  of 
the  medium,  is  under  current  investigation.  This  will  allow  wave  propagation  through 
complex  realistic  media  to  be  stochastically  modeled. 

The  addition  of  these  features  should  elevate  the  phase-screen  method  to 
the  point  where  it  will  be  a  viable  tool  for  seismic  discrimination. 
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